rm(list = ls())
gc()

library(lubridate)
library(tidyverse)
library(foreign)
library(ggplot2)
library(ggthemes)
library(estimatr)
library(PRROC)
library(pROC)
library(glmnet)
library(tmap)    
library(mapview)
library(tmaptools)
library(sf)
library(spData)
library(raster)

source("functions/pca_fun.R")

df = readRDS("data/baseline_clean_mean.rds")

lb_shape = st_read(dsn = "data/lbn_adm_cdr_20200810/lbn_admbnda_adm2_cdr_20200810.shp")

locations = readRDS("data/locations.rds")

locations$location_district[!locations$location_district %in% lb_shape$admin2Name] %>% unique

df$location_district = NULL

df = df %>% left_join(locations)

# PCA ---------------------------------------------------------------------
leb_inputs = c("head_work_legal", "head_work_days_4", "head_work_hours_4", "head_income",
               "hh_inc_source_aid", "aid_atm", "aid_wfp", "hh_leb_2011", "hh_inc_stable", 
               "aid_chg", "own_items_1", "own_items_2", "own_items_3", "own_items_4", "own_items_5",
               "own_items_6", "own_items_7", "own_items_8", "own_items_10", "own_items_11",
               "expenses_ability", "assets_value", "assets_months", "hh_income",
               #
               "ipl_connected", "ipl_outsider", "ipl_arabic", "ipl_literate", "ipl_politics", "ipl_discuss", 
               "ipl_job", "length_stay", "curfews_never", 
               "no_curfews_now", "ind_public_treat", "ind_auth_treat", "not_detained", "no_disc_housing", 
               "ind_mobility", "ind_mobility_hh",
               #
               "head_not_sick", "head_treated", "person_nosick", "person_treat", 
               "ipl_doctor", "disc_healthcare", "no_need_school", 
               "stability_towns", "stability_current", "scl_not_preventive", "access_legal", "own_items_9",
               #
               "registered", "resident",
               #
               "hh_syr_leb", "ipl_lebanese", "ipl_syrian", "leban_relatives", "ipl_meal")

df[,leb_inputs] %>% summary

df$push_factors = PCA_extract(vars = df[, leb_inputs], cor_type = "pear")[[2]] %>% as.vector()
x = PCA_extract(vars = df[, leb_inputs], cor_type = "pear")[[1]]
x$Structure

options(max.print=1000000)

df_plot = df %>% 
  group_by(location_district) %>% 
  summarise(push_factors = stats::weighted.mean(push_factors, weights))


df_plot = df_plot %>% rename(admin2Name = location_district)
lb_shape = lb_shape %>% left_join(df_plot)

lb_shape = lb_shape %>% rename(`Conditions in Lebanon\n(standard deviations)` = push_factors)

map = tm_shape(lb_shape) + tm_polygons(col = "Conditions in Lebanon\n(standard deviations)", palette = "RdBu") + tm_shape(lb_shape) + 
  tm_text("admin2Name", auto.placement = F, size = .4,
          xmod = c(0, #1
                   0, #2
                   0, #3
                   -1, #4
                   0, #5
                   0, #6
                   0, #7
                   0, #8
                   0, #9
                   0, #10
                   0, #11
                   0, #12
                   0, #13
                   -.6, #14
                   0, #15
                   0, #16
                   0, #17
                   0, #18
                   0, #19
                   0, #20
                   0, #21
                   0, #22
                   0, #23
                   0, #24
                   0, #25
                   0 #65
                   ),
          ymod = c(0, #1
                   0, #2
                   0, #3
                   .5, #4
                   -.4, #5
                   0, #6
                   0, #7
                   0, #8
                   0, #9
                   0, #10
                   0, #11
                   0, #12
                   0, #13
                   .4, #14
                   0, #15
                   0, #16
                   .2, #17
                   0, #18
                   0, #19
                   .7, #20
                   0, #21
                   0, #22
                   0, #23
                   0, #24
                   0, #25
                   0 #65
          )
          ) + 
  tm_layout(legend.title.size = .7,
            legend.text.size = 0.5,
            legend.bg.alpha = 1)

tmap_save(map, "figures/appendix/map_push_factors.pdf", width = 6, height = 5)
tmap_save(map, "figures/appendix/map_push_factors.png", width = 6, height = 5)
